A power-law distribution of phase-locking intervals does not imply critical interaction 
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Neural synchronisation plays a critical role in information processing, storage and 
transmission. Characterising the pattern of synchronisation is therefore of great in- 
terest. It has recently been suggested that the brain displays broadband criticality 
based on two measures of synchronisation — phase locking intervals and global lability 
of synchronisation — showing power law statistics at the critical threshold in a classical 
model of synchronisation. In this paper, we provide evidence that, within the limits 
of the model selection approach used to ascertain the presence of power law statistics, 
the pooling of pairwise phase-locking intervals from a non-critically interacting system 
can produce a distribution that is similarly assessed as being power law. In contrast, 
the global lability of synchronisation measure is shown to better discriminate critical 
from non critical interaction. 



I. INTRODUCTION 

The notion of criticality has been hotly discussed in re- 
lation to its presence in the human brain 1-5]. Support 
for the concept of a critical brain has emerged from com- 
paring brain dynamics at various scales with the dynam- 
ics of physical systems at criticality. Much impetus for 
this line of work has come from the observation of power 
laws, a necessary but insufficient condition for critical- 
ity,in distributions associated with neuronal avalanches 
[6|,|7|, but further evidence has come from the application 
of methods from statistical physics for identifying spatio- 
temporal scaling functions in fMRI [1, long-range 
temporal correlations in amplitude fluctuations of band- 
pass filtered electro/magneto-encephalogram (M/EEG) 
[Tol [Tl1 | as well as universal scaling functions in the ac- 
tivity of individual neurons [l2l . Il3| . Functionally, it 
has been difficult to attribute relevance to these find- 
ings other than by making observations of difference in 
some scaling parameter between different human subject 
populations or with the subject's age. It would therefore 
be of great interest to find evidence of criticality in the 
synchronisation of activity between different brain areas 
i.e. a parameter that has been directly linked with infor- 
mation processing, storage, and transmission [Til. ITsT]. 

A system at, or close to, a critical phase transition has 
been associated with the possibility of rapid reconfigura- 
tions in response to external stimuli 0, [HI • Kitzbichler 
et al. pj], 13 argue that rapid state changes are cru- 
cial for the brain to deal with the environment it meets. 
They suggest that in some situations, an extensive cogni- 
tive effort is required and information transfer needs to be 
maximised between brain regions, and at other, relatively 
quiescent periods, the greater concern is minimising neu- 
ronal wiring costs [HJ . A brain at criticality might allow 
the necessary rapid transitions in functional connectivity 
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to occur quickly [l9j . Werner [jj| indicates that a neu- 
rophysiological system in a critical state is best able to 
learn and remember complex logical rules, by adapting 
its synaptic weights quickly. Meisel et al. [20] suggest 
that local events can spread rapidly through a system in 
such a state, and that remaining at criticality prevents 
the spread both from becoming uncontrollably large, or 
from dying away without effect. A single element hence 
has the ability to affect the entire system, which may be 
crucial to processing external stimuli efficiently [2lj |. 

To assess criticality of synchronisation, Kitzbichler et 
al. [171 ] proposed two measures characterising the pattern 
of synchronisation in a complex system. The first mea- 
sure is the frequency density of phase locking intervals 
(PLI) , which are defined as the periods of time for which 
two oscillators differ in their phase by less than a value 
of 7r/4 in modulus. The phase, here, describes where an 
oscillator is in its cycle, relative to the origin. It evolves 
in the interval [— 7r, n] as the oscillator completes an os- 
cillation. The second measure is the frequency density 
of the change in number of phase locked pairs between 
successive time points (global lability of synchronisation 
or GLS). Both measures are derived from a thresholded 
wavelet-transformed instantaneous phase difference (fur- 
ther introduced in Sections III. 51 and III.6[) . Kitzbichler 
et al. validated the PLI and GLS results by showing 
that in two known models of critical interaction, namely, 
the Ising model [H [H[ and the Kuramoto Model (24 - 
l26j (further discussed in Section III.1|) . these measures 
display power law distributions at the critical threshold 
but not in a decoupled system jl7| . The presence of this 
power law in the PLI and GLS was determined using a 
model selection approach [13, [28| whereby both the power 
law and alternative models (log-normal and exponential) 
are fitted and the best model is decided on the basis of 
the Akaike Information Criterion (formally introduced in 
Section 111.71) . 



Whilst it is true that power law statistics of some ob- 
servable of the system should be evident in a system 
at criticality @, H^-Ull, the point has been made that 
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power laws could result from the superposition of mul- 
tiple pro cesses each with their own characteristic time 
scale |32j | or from the use of thresholds [33[ ■ Given this, 
we ask whether power law distributions in the PLI and 
GLS measures introduced in (l7| are uniquely indicative 
of a system in a critical state. Our approach is to pool 
the phase locked intervals (respectively, the number of 
phase locked pairs between successive time points) of a 
non-critically interacting system of Kuramoto oscillators 
and compare the resulting distributions with those de- 
rived from a critically-coupled system. If this pooling 
produces distributions that, within the limits of a model 
selection approach, cannot be distinguished from those 
of a critically-coupled system then we suggest that this 
approach to inferring criticality is suspect. To do so, we 
consider a system formed from a collection of indepen- 
dent paired oscillators, which we refer to as the Inde- 
pendent Pairs model. The two oscillators making up a 
pair are coupled, having phases evolving according to the 
Kuramoto differential equations (formally introduced in 
Section III. II) . but there is no connection between pairs. 
Each pair can snap into synchronisation at a coupling 
value unique to itself, however, there is no collective or- 
der parameter to unite their progressive synchronisation, 
i.e., this system can have no critical coupling value. 

The paper is organised as follows. After a brief re- 
view of the Kuramoto oscillators (Section Bl.l|) . we derive 
analytically the phase difference between two sine-phase 
coupled oscillators, which makes it possible to generate a 
large number of Independent Pairs, with natural frequen- 
cies drawn from a normal distribution and pair-wise cou- 
pling a free parameter (Section III.2p . After summarising 
the methodology of Kitzbichler et al. (Sections III.3flII.7[) . 
we compare its application to both the Kuramoto model 
and our Independent Pairs model (Sections IIII.2IIIII.3l) , 
revealing the coupling parameters under which PLIs and 
GLSs may give rise to power laws within a model selec- 
tion approach. 



II. METHODS AND MATERIALS 
II. 1. The Kuramoto model 



i in such a system is given by [24Tj26j : 

Oi = Ui + ^SjLiBin(^ - 0$ (1) 

Kuramoto [24[ showed that the evolution of any phase 
6i can be re-expressed using two mean field parameters, 
which result from the combined effect of all oscillators in 
the system. Namely, we may say: 

§i =uJi + Krsm(i/j - Z ) (2) 

where ip is the mean phase of the oscillators, and r is 
their phase coherence, so that: 

re* = l£v* (3) 

3=1 

This crucially indicates that each oscillator is coupled to 
the others through its relationship with mean field pa- 
rameters r and if>, so that no single oscillator, or oscilla- 
tor pair drives the process on their own. The oscillators 
synchronise at a phase equal to the mean field ip, and r 
describes the strength of synchronisation, sometimes re- 
ferred to as the extent of order in the system [4(| l4lj . 
When r = 0, no oscillators are synchronised with each 
other. When r = 1, all oscillators are entrained with 
each other. 

It is easy to see that one solution to Equation [2] is 
r = for all time and coupling, leaving each oscillator to 
evolve independently at its own natural frequency. Using 
a limit of N — > oo, some further deductions can be made, 
including the fact that when the natural frequency distri- 
bution g(tu) is unimodal and symmetric, another solution 
can be found for Oi, with r not equivalent to [24]]. A 
critical bifurcation occurs for sufficiently high coupling, 
resembling a second-order phase transition [42j in which 
the order parameter (here, r) leaves zero and grows con- 
tinuously with coupling (4(| H^]. The coupling at the 
bifurcation is referred to as the critical coupling K c (43[ . 
While the above definition holds for a system of infinite 
size, for a finite system such as that considered in this 
paper, the critical coupling can only be approximated by 
this theoretical value. In Section Hl.41 we will provide an 
operational definition of critical coupling in a finite size 
system. 



The Kuramoto model is a classical model of synchro- 
nisation [H, [35[. It has been widely used to study the 
oscillatory behaviour of biological systems such as the 
sleep and body temperature cycles in humans [3(1 
heart pacemaker cell firing 34l . l36l |37|. n euronal firing 
[H IM [11 and fire-fly flashing @ [MM HI- 

The Kuramoto model describes the phase behaviour of 
a system of mutually coupled oscillators with a set of dif- 
ferential equations. Each of N oscillators in the system 
rotates at its own natural frequency {uii, i = 1, N}, 
drawn from some distribution g(u>). However, it is at- 
tracted out of this cycle through coupling K, which is 
globally applied to the system. The differential equation 
to describe the time evolution of the phase Oi of oscillator 



II. 2. Analytic Phase Difference for the 
Independent Pairs Model 

An independent pair is defined as two coupled oscil- 
lators i and j whose phases evolve according to Equa- 
tion ([T|), namely: 

K 



% - 0j = fa + — (sm(9j - Oi) - sin(0j - 0j)) 



= (wi — ujj) — K (sin(0j -0j)) 
Letting A t3 — Oi — 0j yields: 

A I; = (u)i - Wj) - Ksxn(Aij) 



(4) 
(5) 
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This equation has two solutions depending on whether 



K < uji 



or K >| 



If we let C* = 



K 



"J I UI n ^1 >"* ^ !■ 11 wc 1CU w — 

and D is an integrating constant, then the solution for 
K <| uji — ujj | is: 

(t-D) {^-u 3 )yf{l-C 2 ) \ 



■ 2tan 



'1 - C 2 ] tan 



The solution for K >| 



is: 



A, 



2tan 



VC 2 



-i(^-^)V(c 2 -i) 



c 



(6) 



(7) 



with A an integrating constant. A full derivation is pro- 
vided in the Appendix. After deriving this, the authors 
were made aware that the dynamics of a single pair from 
this model has previously been described in [4J| in rela- 
tion to the interaction between a pendulum suspended 
in a viscous fluid inside a rotating container, and used in 
p5j as a basis for constructing a Lyapunov function. 

The time evolution of Ay is dependent on two param- 
eters: the coupling K, and the difference between the 
natural frequencies of rotation, uji — ujj of the two oscil- 
lators. The selection of these two quantities is crucial to 
further analysis and we look at each in turn. 



II. 3. Natural Frequencies 

The natural frequencies of oscillators in the Kuramoto 
system considered in [i~7j were drawn from a normal dis- 
tribution Af (0,1). As any normal distribution may be 
scaled and shifted so that it is equivalent to one with a 
mean of and a standard deviation of 1, we consider 
that our natural frequencies are also distributed with 
bJi ~ Af (0,1) without loss of generality. If both natu- 
ral frequencies u>i and ujj are drawn in this way, then by 
laws of normal distributions, Ui — uij ~ Af (0, 2). As the 
quantity uji — ujj only is of interest to us in order to cal- 
culate Ay (Equations |5] and [7|) , we draw values from a 
distribution of Af (0, 2) for the Independent Pairs Model. 



II. 4. Coupling Parameter 

The critical coupling parameter was calculated analyt- 
ically by Kuramoto under a certain set of assumptions 
[24j . Namely, if the probability distribution of the nat- 
ural frequencies g(oo) is unimodal and symmetric, and 
the number of oscillators is infinite (N oo), then the 
analytic critical coupling parameter K c is: 



K P = 



7^(0) 

And, in the case of g(w) = N (0, 1): 
2^2 



1.596 



(8) 



(9) 



In any feasible realisation of the Kuramoto model, the 
assumption N — > 00 is not realistic. This means that 
the theoretical value of K c ~ 1.596 is not necessarily the 
precise coupling parameter for which the system reaches 
critical behaviour. Kitzbichler and colleagues [13] de- 
scribe two practical measures characterising the onset 
of synchronisation with increasing coupling. The first is 
the change in the 'effective mean-field coupling strength', 
A(Kr). If the value of Kr exceeds the difference be- 
tween the natural frequency and the mean phase uj^ — ijj 
(in modulus) i.e. |cUj — tp\ < Kr, then oscillator i will 
synchronise to the mean field [46| . Thus the value of K 
at which Kr increases maximally is the coupling value 
at which the greatest number of oscillators are drawn 
into the mean field, i.e., a defining feature of the critical 
point in the system. The second measure is the change 
in the time-averaged number of synchronised pairs N$p 
as the coupling increases, AN$p- Again, this describes 
the point at which the greatest change in synchronisation 
occurs, i.e., the critical point. The two measures A(Kr) 
and AN$p peak at the same point. We shall call the 
coupling value at this point the effective critical coupling 
value for our system. 

In contrast, in our Independent Pairs model, there is no 
longer a global critical coupling parameter K c since there 
can be no mean field. From the two distinct analytical 
solutions for Ay (Equations [5] and [7]) we see that each 
pair of oscillators will synchronise independently when 
K exceeds | uji — ujj \ for that pair. Some insight can 
nevertheless be gained by calculating the measures de- 
rived from a standard Kuramoto model, namely, r, N$p, 
A(Kr) and AJV SP . 

As shown by Figure [TJA, there is a clear growth in or- 
der in the Kuramoto model, with the parameter begin- 
ning near for low coupling, and increasing to nearly 1 
after the coupling value exceeds K — 3. The maximum 
rise in Kr occurs at around K = 2, which is therefore 
the effective critical coupling for this system. A similar 
pattern is traced by N$p, with AN$p peaking at around 
K = 2. In this paper, we will provide results for the theo- 
retical critical value K c ~ 1.596 (occasionally referred to 
as K c ~ 1.6), as well as for the (above defined) effective 
critical coupling for our finite system, K = 2. This latter 
value is where we might expect power law statistics to 
be present in the Kuramoto model. The authors have 
empirically confirmed that as N increases, the effective 
critical coupling K converges to the theoretical critical 
coupling K c (results not shown, but the effective critical 
coupling is K = 1.8 for N = 1000 for example). It should 
be noted that although the number of oscillators consid- 
ered here is limited, 44 oscillators as in [l?]], this system 
still gives rise to 946 pairwise interactions, which is more 
substantial. From a neuroscience viewpoint, it could be 
argued that 44 oscillators are sufficient for drawing useful 
conclusions about a neuronal system. For example, the 
use of a Kuramoto model of 66 phase oscillators by the 
authors of [13] led to the emergence of slow activity fluc- 
tuations consistent with empirically measured functional 
neural connectivity. Nevertheless, in order to verify our 
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A. Kuramoto Model 
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B. Independent Pairs Model 
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FIG. 1. Plot A. shows the evolution of order parameter r for the Kuramoto model with cyan solid circles (error bars show 
standard deviations). The coupling parameter K increases along the a;-axis. The hollow purple diamonds show AKr, the 
change in order parameter multiplied by coupling (error bars not shown for readability). The time-averaged number of 
synchronised pairs, Nsp is shown with hollow green squares (error bars show standard deviations), and the difference in Nsp, 
ANsp, is indicated by solid blue triangles (error bars not shown for readability). The peaks in AKr and ANsp can be used to 
indicate the location of the critical point for a specific system, which for this selection of natural frequencies occurs at around 
K = 2. This effective coupling value of K = 2 will be used throughout the paper. Note that, for the Kuramoto model, the 
order increases with rising coupling. Plot B. displays the corresponding measures r, AKr, Nsp and ANsp for the Independent 
Pairs model. There is no change in order parameter with coupling, indicating that the oscillators are not critically coupled to 
a mean field. 



conclusions, we replicated our analysis with N = 1000 os- 
cillators yielding similar results (not shown but available 
upon request from the corresponding author). 

With independent pairs, on the other hand, both the 
order parameter and the number of synchronised pairs 
remain unchanged across all coupling values, at the val- 
ues observed for K = in the Kuramoto model (see 
Figure QJ3). This is because, although the pairs individu- 
ally synchronise with each other, the frequencies at which 
they synchronise are distributed across the whole range 
of possible frequencies. 



II. 5. Frequency scales 



An important feature of the findings in [17j is that 
the critical behaviour of neural activity extends across a 
number of frequency scales, so that criticality is referred 
to as being broadband. The decomposition of the phase 
difference data into several frequency scales is done us- 
ing a Hilbert wavelet transform, and was implemented 
computationally here using the algorithms from (4cj - 
[50j . Specifically, wavelet scales 3-11 were used, cor- 
responding to frequencies of 125 — 62.5Hz, 62.5 — 31Hz, 
31 - 15.5Hz, 15.5 - 8Hz, 8 - 4Hz, 4 - 2Hz, 2 - 1Hz, 
1 - 0.5Hz, and 0.5 - 0.25Hz. 



First, Kitzbichlcr et al. [17[ construct two signals 
denoted Sj and Sj hereafter, by taking the cosine of 
phases 6i and 9j respectively. They then take the k-th 
scale wavelet transforms of s, and Sj to obtain Wfe(si) 
and Wfc(sj), which are time- varying complex vectors 
of wavelet coefficients. Each set of wavelet coefficients 
quantifies the power of the signal in the correspond- 
ing frequency band. These two sets of wavelet coef- 
ficients are multiplied element-wise to form the vector 
Wk{si)^Wk{sj), where the symbol f indicates the com- 
plex conjugate. This vector is then normalised by divid- 
ing it (again, element-wise) by the element-wise product 
I Wfc(sj) || Wfc(sj) | where operator | . | denotes the 
modulus. The result is an instantaneous time-varying 
complex phase vector: 



W k {s^W k {s 3 ) 
Wk(si) || W k (sj) | 



(10) 



To ensure a more robust and less noisy estimate of 
the phase relation, the instantaneous phase vector is 
smoothed by using a moving average of the numerator 
and the two vectors contributing to the denominator of 
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Cfij, yielding a new vector Cy given by: 

13 y/(\WM P)P^IP) 

Here the operator (.) denotes that a moving average is 
taken. The length of the sliding window used for the 
moving average is set to the number of time steps span- 
ning 8 oscillation cycles at the highest frequency in that 
wavelet scale 

The argument of Cy is then taken as a measure of the 
phase relationship of the two oscillators i and j corre- 
sponding to wavelet scale k, so that Ay = arg(Cy). 

In the Independent Pairs model, the phase differences 
within each pair are known analytically (see Section lII.2p . 
however, they are not associated with particular wavelet 
scales. To produce probability distributions comparable 
to those in [l7|, surrogate pairs of signals were created 
with the first signal evolving constantly at a frequency 
given by a base value drawn from the distribution of nat- 
ural frequencies g(ui), and the second signal phase shifted 
from the first by Ay. 

II.6. PLI and GLS 

In this section, we will use Ay (t) to denote the value of 
A^ at time t. For phase difference Ay. between two os- 
cillators i and j, the PLIs are defined as the duration (in 
seconds) for which —a < Ay-(i) < a, for some threshold 
a. This definition was given by 17] with a — ir/4. 

The GLS was also defined in [17[ and characterises the 
evolution of the number of synchronised pairs, Nsp, to 
describe the lability of synchronisation. The number of 
synchronised pairs at wavelet scale k is formally defined 
as: 

N h SP (t) = {l A U*) l< « and M V® > 1} (!2) 

where M^ 2 =| Ch | 2 is proposed as a measure of the 
significance of the phase difference estimate Cy , and a — 
7r/4 as above. It should be noted here that the condition 
My (t) > | introduces an additional threshold. The 
use of thresholds on otherwise stochastic data has been 
shown by Touboul et al. [33[ to occasionally give rise to 
spurious power laws. 

The GLS at scale k is then obtained by calculating the 
square of the difference in the number of phase-locked 
pairs between two successive points in time: 

GLS k =| Ng P (t + St) - N$ P (t) | 2 (13) 

where St is an increment in time and k denotes the 
wavelet scale. 

From examination of our analytic equations for phase 
difference (Equations[6]and[7|), we observe that the phase 
difference Ay changes with time in a very structured 
way. For K < | u>i — ujj |, Ay is a periodic function. For 



K > | LJi — LUj | , there is a short-lived transient before Ay- 
settles to a constant. 

Before we proceed to pool our probability distributions 
across many pairs of oscillators, we first consider what we 
might expect from a single pair. 

For K < | u!i — uij | , the lengths of PLIs between two 
oscillators would be identical within any given oscillation 
cycle, and the probability distribution will only contain 
one value. If a given simulation is cut off before a full 
cycle is complete, or more precisely, before a phase locked 
interval has come to an end, this may give rise to a second 
phase locked interval, and the probability distribution 
may have more than one value in this case. For K > \ LOi — 
Ljj |, the phase difference will be a single constant, either 
occurring during the transient, or at the permanent value 
to which the phase difference converges, depending on 
the starting phase difference, and the value of the final 
constant. Again, the probability distribution contains 
one value. 

The GLS can either take the value 1 if the oscillators 
either go from being non-phase-locked to phase locked, or 
the value if no change occurs. This allows two possible 
values in the probability distribution. 

For a single oscillator pair, we would therefore not ex- 
pect to find a valid probability distribution of either PLIs 
of GLS for any coupling K. 

This is a trivial, but important point to make. If a 
single pair of oscillators could give rise to a probabil- 
ity distribution which appeared linear on a log-log plot 
(as a power law does) for some pairwise coupling value 
that could be considered 'critical' over some small range 
of values, then the final, observed power law created by 
pooling many pairs may be the result of a simple su- 
perimposition of these smaller linear components. We 
now demonstrate that the power law could result from 
a process that does not involve 'critical' interactions for 
any reasonable definition of the term (even on a pair- 
wise level) , but through completely independent systems 
evolving with no connections between the elements that 
combine to produce the power law. 



II. 7. Akaike Information Criterion 

As in [T3|, the presence of power law statistics is as- 
sessed using a model selection approach whereby the 
Akaikc's Information Criterion [5l| is used to compare 
the goodness-of-fit of a power law distribution with that 
of two alternative distributions, namely, the exponential 
and log-normal distributions. It is important to stress 
that the Akaike Information Criterion only provides a 
means of comparing models, but gives no information on 
how good the model is objectively at fitting the data. 
This means that only the relative values of this measure, 
for different models, are important. 

For a model using k parameters, with likelihood func- 
tion L, the Akaike Information Criterion is calculated 
using the following expression: 
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A. C = 0.00 



B. C = 0.50 



AIC = 2k- 21n(L) 

As in this measure was adjusted to account for 
small sample sizes, using the following: 



AICc = AIC 



2k(k + 1) 



n — k — 1 

where n is the number of observations of the data. 
This is especially relevant because all three models were 
fitted to the binned histogram heights, rather than the 
full data set. Since the basis of the AIC is a log-likelihood 
function, it can be used with binned data in this way [52j ■ 
The number of bins used will affect the raw values of the 
AIC, but not the relative values obtained for the models 
used, so that the best-fitting model will pertain for the 
data analysed. 



III. RESULTS 
III.l. Independent Pair model simulation 

We simulated pairs of Kuramoto-coupled oscillators 
alongside our analytic solution. Both were calculated 
over 1000 seconds, with an integration time step of 
St = 2~ n for the simulated oscillators. This provided 
a total of 1000 x 2 11 time steps. We then down-sampled 
the resulting time series by a factor of 2 to obtain a time 
series with sampling frequency of 2 10 Hz. The analytic 
signal was also generated with a sampling frequency of 
2 10 Hz. The coupling K was incremented between and 
4, in intervals of 0.2, and the two curves were compared. 

The behaviour of the phase difference is qualitatively 
different in the cases C = 7 K . < 1 and C > 1. 
We demonstrate the phase difference between two os- 
cillators in Figure [5] as obtained with our analytic ex- 
pressions alongside a simulation of the Kuramoto model, 
using Euler's method to iteratively update the phase by 
Equation [TJ The two phase calculations are perfectly su- 
perimposed. 

Although the root mean square error (RMSE) varies 
for different coupling values, the normalised RMSE is less 
than 0.1% for the range of coupling values considered 
in this paper, demonstrating good agreement between 
simulated and analytic results. 

It is evident that when the coupling supersedes the 
difference in natural frequencies (C > 1), the two oscilla- 
tors synchronise in exponential time. When the coupling 
is small (C < 1), however, the phase difference grows 
(or falls) at a rate dictated by the frequency difference, 
but with increasingly lengthy periods of constant phase 
difference, or synchronisation. 



III. 2. PLI and GLS of Kuramoto model 

As a baseline for comparison, the results of Kitzbich- 
ler et al. [Tz| on the Kuramoto model were replicated 
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FIG. 2. The evolution of phase difference between the oscil- 
lators in a two-oscillator Kuramoto system, plotted using our 
analytic expression (blue), and a simulation of the Kuramoto 
model by Euler's method (red). Fhe two phase calculations 
are perfectly superimposed. The root mean square error 
(RMSE) is shown for different coupling values, for a single 
simulation. Panels A,B,C have C < 1 (where C is defined 
in Section [II. 2\ , but coupling is increased progressively. The 
phase evolves periodically. Panel D is the same pair of 
oscillators, but for C > 1. There is a brief transient before 
the oscillators fully synchronise with a constant level of 
phase difference. The initial phase separation has been set to 
A = without loss of generality. 



using our own code in the Matlab environment. A sys- 
tem of 44 Kuramoto oscillators, each with a natural fre- 
quency drawn from a normal distribution M (607T, 207r), 
was simulated using the same simulation parameters as 
in Section lilI.il We present three different regimes (un- 
coupled, critically coupled, and super-critically coupled), 
which yield the power spectra shown in Figure [3J 

Next, using 44 oscillators whose natural frequencies 
were drawn from a M (0,1) distribution, the PLI and 
GLS probability distributions were calculated for the fol- 
lowing coupling values - K = 0, K = K c = 1.596, K = 2 
and K = 4. At t = 0, all oscillators had a phase ftj = 0. 
The data presented in figures |3l [4] and [5] were obtained 
from a single run of the model, however, it was confirmed 
that the results were not sensitive to the exact values of 
the natural frequencies. 

A histogram for the PLI data was constructed using 20 
logarithmically spaced bins, with the first bin beginning 
at a single time step of 2~ 10 seconds, and the largest bin 
ending at the total length of the data, of 1000 seconds. 
The histogram was then scaled so that each bin count 
was divided by the total number of PLIs, and then by 
the bin size that it represented. 
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Coupling: 0.00 Coupling: 2.00 Coupling: 4.00 




frequency (Hz) frequency (Hz) frequency (Hz) 

FIG. 3. Power spectra for a system of 44 Kuramoto oscillators, with natural frequencies drawn from a M (607r, 207r) distribution 
and three distinct levels of coupling - A) K = 0, B) K = 2, the effective critical coupling for this specific finite Kuramoto 
system, as seen from Figure [TJ\ and C) K = 4. The vertical numbered lines represent wavelet scales 3 — 11. 



For GLS, we took 1000 logarithmically spaced bins 
ranging from a value of 1 to 10 4 5 , as displayed on the 
plot. The GLS histogram was also scaled. Here each bin 
count was divided by the total number of counts (sum 
of all bin counts), and then by the bin size that it repre- 
sented. 

TABLE I. Akaike Information Criterion values for various 
models applied to the PLI distributions of the Kuramoto 
model at K = 2, the effective critical coupling value for our 
system. Smaller values indicate a better fit, but comparisons 
are only meaningful across rows. The smallest value in each 
row is indicated with an asterisk. 

Wavelet Scale Power-Law Exponential Log-Normal 



3 


251.04 




288.75 


116.26 


* 


4 


253.87 




289.35 


123.10 


* 


5 


257.03 




316.55 


157.24 


* 


6 


258.62 




370.14 


218.44 


* 


7 


254.59 




396.20 


252.47 


* 


8 


245.74 


* 


359.41 


250.97 




9 


220.50 


* 


343.30 


227.93 




10 


224.56 




318.80 


229.26 




11 


220.38 


* 


306.27 


223.93 





TABLE II. Akaike Information Criterion values for various 
models applied to the GLS distributions of the Kuramoto 
model at K = 2, the effective critical coupling value for our 
system. Smaller values indicate a better fit, but comparisons 
are only meaningful across rows. The smallest value in each 
row is indicated with an asterisk. 



Wavelet Scale Power-Law 




Exponential 


Log-Normal 


3 


-2533.43 


* 


-1019.49 


-2478.83 


4 


-2531.41 


* 


-1296.02 


-2484.28 


5 


-2540.75 


* 


-1351.52 


-2490.46 


6 


-2520.30 


* 


-1304.60 


-2473.17 


7 


-2439.44 




-1293.77 


-2465.53 * 


8 


-2415.82 




-1163.59 


-2426.63 * 


9 


-2000.55 


* 


-941.78 


-1985.62 


10 


-1536.79 


* 


-686.48 


-1515.75 


11 


-546.67 




-239.38 


-568.82 * 



The Akaike Information Criterion (AIC) was calcu- 
lated for both the PLI and GLS distributions for all stud- 
ied coupling values. Only PLI intervals of length 0.1 sec- 
onds or more were used for model-fitting, and these only 
are shown in the plot. The power-law model was fit- 
ted using the procedure described by Clauset et al. (H^ . 
and implemented using their freely available code, and a 
minimum data value of 0.1 seconds. The log- normal and 
exponential distributions were both fitted using in-built 
Matlab functions. 

The values obtained for the effective critical coupling 
K = 2 are shown in TableUfor PLIs and Table|n]for GLS. 
As in [17], the power law distribution was only found to 



be the best fit at certain wavelet scales. The AIC values 
in Table 1 of Kitzbichler et al. [l7|, stated as being at 
critically coupled Kuramoto, favour a power law model of 
the PLI frequency distribution for 5 of 9 wavelet scales, 
although no value is reported for wavelet scale 11. 
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Coupling: 0.00 Coupling: 1 .60 

10 B n 




PLKs) PLKs) 

FIG. 4. Distribution of PLIs in a system of 44 Kuramoto 
oscillators, with natural frequencies drawn from a M (0, 1) dis- 
tribution and four levels of coupling - K — 0, K — K c ~ 1.6, 
K — 2 and if = 4 (from top-left, clock- wise). A power law 
of exponent -2 is shown by a dotted black line. The coloured 
lines represent wavelet scales 3 — 11 (see key). 



Coupling: 0.00 Coupling: 1 .60 




Coupling: 2.00 Coupling: 4.00 




GLS GLS 

FIG. 5. Distribution of GLS in a system of 44 Ku- 
ramoto oscillators, with natural frequencies drawn from a 
N (0, 1) distribution and four levels of coupling - K = 0, 
K — K c ~ 1.6, K — 2 and K = 4 (from top-left, clock-wise). 
A power law of exponent -1 is shown by a dotted black line. 
The coloured lines represent wavelet scales 3 — 11 (see key). 

In our system, at the effective critical coupling K = 2, 
the power law distribution was the best model for the 



data for 4 out of 9 wavelet scales for the PLI data. Note 
that the same number of wavelet scales were also best 
fitted by a power law distribution for coupling values 
K = 1, K = 3 and K = 4. At coupling K = K c = 1.596, 
3 wavelet scales were best fitted by a power law, and at 
no coupling, i.e., K = 0, only 2 wavelet scales. The log- 
normal distribution was otherwise the best fit at all cou- 
pling values and all other scales. The fact that less than 
half of the wavelet scales were best fitted by a power law 
distribution at the critical coupling, combined with the 
fact that non-critical coupling parameters {K = 1,3,4) 
resulted in the same proportion of scales being best fitted 
by a power law distribution, leads us to conclude that the 
distribution of PLIs is not a reliable measure of criticality 
in a finite size Kuramoto system. 

For the GLS probability distribution the coupling val- 
ues giving greatest resemblance to power law distribu- 
tions were K — K c ~ 1.6 and also K — 3, both with 8 of 
9 wavelet scales best fitted by the power law model. (The 
AIC values for the GLS distribution were not included 
in [13]). In contrast, a power law model was best-fitting 
for only 2 wavelet scales at coupling value of K = 0. It 
was the best fit for 4 wavelet scales at coupling K = 1, 
for 6 wavelet scales at coupling K — 2 and for 3 wavelet 
scales at coupling K = 4. The remaining wavelet scales 
for all coupling values were again best fitted by a log- 
normal distribution. The prevalence of good power law 
fits in the GLS probability distribution across wavelet 
scales for coupling values K = K Cl 2 and 3, and the fact 
that power law distributions were not a good fit for the 
data resulting from coupling values K = and K = 4, 
collectively suggest that the GLS measure may be an ac- 
ceptable but not very sensitive indicator of the region of 
critical coupling for the finite size Kuramoto system. 

The probability distributions of PLIs and GLS in Fig- 
ures EO and [5] are consistent with those shown in Figure 
3 of [17| for the zero and critical coupling values. For 
K = 0, the probability distribution of the PLIs has a 
drop-off for PLI values above 10°. However, our plo t 
at this value differs from that in Kitzbichler et al. [17j . 
which shows that no intermediate length PLIs exist for 
many of the scales. We observe PLIs of all lengths from 
0.1 to over 100 seconds with non-zero probability. We 
suspect that their data was truncated for display, but 
no detail is given in the paper. The distributions at all 
wavelet scales appear linear in the log-log space both at 
theoretical critical coupling of K c ~ 1.6, and at K = 2, 
the effective coupling parameter for this simulation of 
the Kuramoto system. The range in which this linearity 
holds is similar to that in [l7| . lying between 10° and 
10 2 . Our results for coupling values beyond criticality 
show that the distributions remain power-law-like as the 
coupling is increased to K = 3, suggesting that linear- 
ity in the log-log space is not specific to K — K c for 
this system. This linearity in the log-log space vanishes 
for K — 4, where sufficiently many oscillators have syn- 
chronised at the mean field phase for the system, which 
induces a particular interval of phase-locking, indicated 
by the peak in the distribution. Qualitatively similar ob- 
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servations can be made regarding the GLS distributions. 



III. 3. PLI and GLS in the Independent Pairs model 

PLI and GLS probability distributions were computed 
from the phase difference of 1000 pairs of oscillators with 
u>i — uJk ~ Af(0, 2). The length of data, and time steps 
used were identical to those described in Section IIII.ll 
The number of pairs was set to a value close to that of 
the total number (946) of pairings available in a system 
of 44 oscillators. We computed all PLIs across these pair- 
ings, and the measures of GLS for all consecutive time 
points. Histograms of PLI and GLS, and AIC values 
were computed exactly as in the previous Section (see 
Figures © and and Tables El and EJ). 

III. 3.1. PLI probability distribution 

As indicated by Figure HI the structure of the proba- 
bility distribution alters as the coupling increases. For 
K = 0, there is a drop-off below the power law of the 
distribution for values of the PLI above 1 second. At 
or around the theoretical and effective critical couplings, 
the log-log plot of the distribution approaches the same 
power law with slope —2 as indicated by [17]. For values 
up to K = 3, there is no significant difference between the 
evolution of PLI probability distributions with coupling 
in the Independent Pairs model and that of the Kuramoto 
model. The main dissimilarity arises from the continu- 
ing presence of an apparent power law distribution in the 
'super-critical' range of K = 4. In the Independent Pairs 
model, the log-log plot of the distribution retains some 
of its linearity whereas there is synchronisation to the 
mean field in the Kuramoto model, as evidenced by a 
well-defined peak in Figure H) 

For the Independent Pairs Model, the AIC indicated 
that the power law distribution best fitted the PLI proba- 
bility distribution for 4 of the 9 wavelet scales, at critical 
coupling value K ~ 1.6, as well as for coupling values 
K = 1 and K = 4. Both the effective critical cou- 
pling value K — 2 (see Table llPI) and K = 3 favoured 
the power distribution for 5 wavelet scales in contrast to 
only 1 wavelet scale for coupling K = 0. The remain- 
ing wavelet scales at all coupling values were best fitted 
by a log-normal distribution. As there is little differ- 
ence between the numbers of wavelet scales best fitted 
by a power law distribution for corresponding coupling 
values of the Kuramoto and Independent Pairs models, 
we conclude that the PLI measure is therefore unable to 
distinguish between critically and non-critically coupled 
systems. 



III. 3. 2. GLS probability distribution 

In contrast to the PLI results, the probability distribu- 
tion for the GLS of the Independent Pairs model remains 



Coupling: 0.00 Coupling: 1 .60 

10 6 n . 




10 I I 

0.01 0.1 1 10 100 1000 0.01 0.1 1 10 100 1000 

PU(s] PU(s) 



FIG. 6. Distribution of PLIs in the Independent Pairs Model, 
with natural frequencies drawn from a A/"(0, 1) distribution 
and four levels of coupling - K = 0, K = K c ~ 1.6, K — 2 
and K = 4 (from top- left, clock- wise). A power law of 
exponent -2 is shown by a dotted black line. The coloured 
lines represent wavelet scales 3 — 11 (see key). 



Coupling: 0.00 Coupling: 1.60 




Coupling: 2.00 Coupling: 4.00 




GLS GLS 

FIG. 7. Distribution of GLS in the Independent Pairs Model, 
with natural frequencies drawn from a A/"(0, 1) distribution 
and four levels of coupling - K = 0, K = K c ~ 1.6, K — 2 
and K = 4 (from top- left, clock- wise). A power law of 
exponent -1 is shown by a dotted black line. The coloured 
lines represent wavelet scales 3 — 11 (see key). 
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TABLE III. Akaike Information Criterion values for various 
models applied to the PLI distributions of the Independent 
Pairs Model at K — 2, the effective critical coupling value for 
our system. Smaller values indicate a better fit, but compar- 
isons are only meaningful across rows. The smallest value in 
each row is indicated with an asterisk. 

Wavelet Scale Power-Law Exponential Log-Normal 



3 


205.74 




121.02 


49.49 


* 


4 


189.05 




222.37 


120.70 


* 


5 


171.14 




192.08 


107.80 


* 


6 


154.09 




166.67 


93.89 


* 


7 


138.37 


* 


241.74 


139.03 




8 


122.33 


* 


210.90 


124.66 




9 


104.09 


* 


174.94 


109.51 




10 


88.21 


* 


161.30 


93.26 




11 


72.94 


* 


129.74 


80.59 





largely unaltered as coupling increases, as shown in Fig- 
ure [7] The GLS distributions do not resemble those of 
the Kuramoto model. The range in which the log-log plot 
of the distribution is linear is narrower with a drop-off in 
the distribution for values of GLS above 100s, suggesting 
that the Global Lability of Synchronisation measure may 
be more sensitive to the lack of critical interaction in the 
system. 

For GLS, only 2 wavelet scales were best modelled by 
the power law model at the effective critical coupling K = 
2 (see Table HVl for K = K c ). 1 wavelet scale was best 
fitted by a power law at coupling K = 0, 3 at K = 
1, 2 at K = K c , 4 at K = 3, and 3 at K = 4. The 
remaining wavelet scales at all coupling values were best 
fitted by a log-normal distribution. There is no evident 
pattern of increasing similarity to a power law of the GLS 
distribution, as the coupling increases. 



TABLE IV. Akaike Information Criterion values for various 
models applied to the GLS distributions of the Independent 
Pairs model at K = 2, the effective critical coupling value for 
our system. Smaller values indicate a better fit, but compar- 
isons are only meaningful across rows. The smallest value in 
each row is indicated with an asterisk. 

Wavelet Scale Power-Law Exponential Log-Normal 



3 


-297.16 


42.78 


-301.51 


* 


4 


-379.92 


8.93 


-391.39 


* 


5 


-591.87 


-54.62 


-596.56 


* 


6 


-409.53 


-38.71 


-425.36 


* 


7 


-227.94 


-6.39 


-251.63 


* 


8 


-193.42 


23.66 


-204.54 


* 


9 


-129.49 


51.58 


-132.82 


* 


10 


-84.46 * 


57.75 


-78.53 




11 


-63.34 * 


62.20 


-51.41 





IV. CONCLUSIONS 

In this paper, we critically examined two measures, 
phase-locking intervals (PLI) and global lability of syn- 
chronisation (GLS), proposed by Kitzbichler and col- 
leagues [13] to characterise the presence of critical syn- 
chronisation in a system. We did so by presenting those 
measures with two very different models of synchronisa- 
tion. In the first (Kuramoto Model) the oscillators are 
coupled with increasing K to the mean field and undergo 
a critical transition. In the second (Independent Pairs 
Model) the oscillators are only allowed to couple in a 
pair wise manner. This latter model cannot be formu- 
lated as a system at criticality because there is no global 
coupling to associate the pairs with one another, and so 
no possibility of a mean field. 

When calculating the phase locking intervals (PLI) fol- 
lowing the methodology of Kitzbichler et al. [l7|, we 
showed that power laws were the best fit for a similar 
number of wavelet scales when considering PLI distri- 
butions for the critical, Kuramoto, model and the non- 
critical, Independent Pairs, model. The power law dis- 
tribution and the slope found for the PLIs of the non- 
critical system was closely similar to that shown by the 
critical model. When further exploring the PLI probabil- 
ity distribution for coupling parameter values exceeding 
criticality, we found that the linearity of the log-log plot 
of the distribution at a number of wavelet scales still led 
to a best fit by a power law, suggesting that the observa- 
tion of power laws within this framework can be present 
in a wide range of coupling values. We therefore con- 
clude that the PLI measure should not be used to infer 
criticality (broadband or otherwise) in a system. 

In our simulations the GLS measure appeared better 
at discriminating between the critical, Kuramoto, sys- 
tem and the non-critical, Independent Pairs, model. We 
therefore conclude that GLS is a better measure than 
PLI for identifying critical systems, however, we believe 
that further work should be done to ascertain more pre- 
cisely where its strengths lie, and compare it to other, 
non threshold-based methods such as proposed by Gong 
et al. [54]. In particular, we note that the GLS measure 
relies on counting the number of synchronised oscilla- 
tors and that this depends crucially on how oscillators 
are defined, and distinguished. In the Kuramoto model, 
the number of oscillators is well defined, and each one 
is a discrete entity. With recorded neural activity, how- 
ever, distinguishing multiple discrete oscillators is less 
straightforward. Kitzbichler et ah, have applied the GLS 
measure to fMRI and MEG signals but its interpretation 
was limited by finite size effects (see loss of log-log linear- 
ity in the GLS distribution of MEG data in their figures 
5D and 7D). To our knowledge the GLS measure has 
not been applied again to human neural data. Recently 
Meisel et al. have claimed to detect when compared 
to seizure-free electro-corticogram (ECoG) data a loss 
of adaptive self-organized criticality of the ECoG during 
epileptic seizures. This conclusion was arrived at through 
exploring power law scaling of ECoG phase locking using 
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the PLI measure only. This is an exciting finding which 
received support from analysing the changes in PLI scal- 
ing seen in a computational model of self-organized criti- 
cality [55|. However, our work indicates that interpreting 
the presence of a power law in the PLI probability distri- 
bution as a marker of criticality is problematic especially 
when a threshold has been applied to detect PLIs and 
when there has been pooling across many elements. 
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We can solve this integral using the fact that 
tan -1 (z) = J to get: 



fa-wjWO. -CP) 



tan 



tan 



-c 



—tan 



tan 



-c 



VI -C 2 



(A2) 

Here, A?- is the value of Ay at time t = 0, i.e., the 
initial difference in phase between oscillators i and j. 

tan( ^ji )-c 



Setting D 



=tan 



we can rearrange Equation IA2I to get: 




Appendix A: Analytic Derivation of Ay 

The analytical solutions for A y -, the difference between 
phases 0i and 9j of oscillators i and j , are distinct for the 



two cases 



K 



> 1 and 



A" 



< 1 where u>i and uij are 



the respective natural frequencies of oscillators i and j, 
and K is the coupling added globally to the system. We 
can rearrange Equation[5]to obtain the following integral: 



dt 



dA 



(iJi - Ljj) - Ksin(Aij) 



where t denotes time. This integral can be solved using 
the standard substitution of x = tan {^^j ■ 



Doing so, and letting C — ( u f_ u .y we get: 



dt = 



dx 



(ui-Uj)J (1 - CP + (x - C) 2 ) 



(Al) 



There are two different scenarios for this integral, de- 
pending on whether C < 1 and vl — C 2 is a real or 
imaginary number. We deal with each case in turn. 



2. If C > 1, or when coupling is larger than the 
difference in natural frequency 



Here, vT — CP is imaginary, so we rearrange IA1I in 
terms of \J C 2 — 1: 



dt 



( Wl -^)(1-C 2 ) 



dx 



1 - 



x-C 



We can solve this integral using the fact that 
\ (log-^-z - 1) - log-^z - 1)) = / j^: 



t 



-1 



(uJi-cjjWiC 2 -!) 



log 



.4 



l+y 



where A — j+p- and y° is the value of y at time t = 0. 
This can be rearranged to yield: 



Aij = 2tan _1 



, / p -tK- W3 )V((7 2 -i) _ A s 

VC^ll- ; ~ 



c 



1. If C < 1, or when coupling is smaller than the 
difference in natural frequency 



We can rearrange I All in terms of a/1 — C 2 which is real 
and: 
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